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ABSTRACT 

Because of their simplicity, axisymmetric mass distributions are often used to model 
gravitational lenses. Since galaxies are usually observed to have elliptical light dis- 
tributions, mass distributions with elliptical density contours offer more general and 
realistic lens models. They are difficult to use, however, since previous studies have 
shown that the deflection angle (and magnification) in this case can only be obtained 
by rather expensive numerical integrations. We present a family of lens models for 
which the deflection can be calculated to high relative accuracy (10 -5 ) with a greatly 
reduced numerical effort, for small and large ellipticity alike. This makes it easier to use 
these distributions for modelling individual lenses as well as for applications requiring 
larger computing times, such as statistical lensing studies. A program implementing 
this method can be obtained from the author^. 

Subject headings: gravitational lensing — cosmology — galaxies: structure 



1. Introduction 

The types of density profiles used to model gravitational lenses have been motivated by obser- 
vations of lenses in addition to practical considerations. Observed features of galaxies and clusters 
that can be incorporated into lens models include ellipticity and radially decreasing mass density 
profiles. But it is essential to be able numerically to calculate quickly the deflection angle and 
magnification of a light ray due to a lens model, in order to probe the entire range of parameter 
space when searching for a best fit to an observed lens. Numerical efficiency is even more important 
in cases with multiple sources such as radio lenses observed at high resolution. Strongly lensed arcs 
in clusters sometimes lie near cluster galaxies and thus modelling all the lensed features in a clus- 
ter may require a combination of many individual lenses. Another application which depends on 
numerical speed is the statistical study of properties of the images of a lens model for comparison 
with lens surveys. 
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In attempting to construct models that are as realistic as possible, new features must be 
introduced carefully, since in a single case of multiple lensing there are a small number of constraints 
requiring a small number of model parameters. Ellipticity adds just two parameters (magnitude 
and orientation), and it appears to be essential. Galaxies and clusters often appear elliptical, 
with significant numbers having axis ratios b/a smaller than 0.5. Axisymmetric lens models are 
also excluded by many observed lens systems, in which the images are not colinear with the lens 
position on the sky. Of course, the asymmetry in each case may also be due in part to external 
shear from other nearby galaxies or from large-scale structure along the line of sight (Bar-Kana 
1996, Keeton et al. 1997). 

However, the use of elliptical mass distributions has not become common practice because the 
evaluation of the deflection angle and magnification matrix requires some effort. Bourassa et al. 
(1973) and Bourassa and Kantowski (1975) (with minor corrections by Bray (1984)) introduced a 
complex formulation of lensing which allows for an elegant expression of the deflection angle due to a 
homoeoidal elliptical mass distribution. I.e., this is a projected, two dimensional mass distribution 
whose isodensity contours are concentric ellipses of constant ellipticity and orientation. It can 
be obtained, e.g., by projecting a three-dimensional homoeoidal mass distribution. The complex 
integral which gives the deflection angle is in practice difficult to separate into real and imaginary 
parts. Schramm (1990) used an alternative derivation to obtain the deflection without the use 
of complex numbers, but still requiring a numerical integral for each component of the deflection 
angle. Elliptical densities have been used for numerical lens modelling, with the ellipticity allowed 
to vary in order to fit the data, by Keeton & Kochanek (1997). 

In order to avoid numerical integration, several alternatives have been suggested to exact el- 
liptical mass distributions. Models where the potential is chosen to have elliptical contours rather 
than the density are easy to use, since the deflection can be obtained immediately as the gradient 
of the potential. The imaging properties of elliptical potentials have been investigated extensively 
(Kovner 1987, Blandford & Kochanek 1987 and Kochanek & Blandford 1987). They become iden- 
tical to elliptical densities for very small ellpticities and produce similar image configurations even 
for moderate ellipticty (Kassiola & Kovner 1993, hereafter KK93). However, elliptical potentials 
cannot represent mass distributions with axis ratios b/a smaller than about 0.5 because the cor- 
responding density contours acquire the artifical feature of a dumbbell shape, and the density can 
also become negative in some cases (Kochanek & Blandford 1987, KK93). To avoid this problem, 
Schneider & Weiss (1991) proposed a numerical method based on a multipole expansion of the 
mass distribution, but this expansion converges slowly when used with large ellipticity. 

In this paper we consider a family of projected density profiles that has been used with the 
approximate approaches to ellipticity discussed above. This is the family of softened power-law 
profiles, which have a constant density within a core radius and approach a power-law fall-off at 
large radii. In §2 we introduce our notation for softened power-law elliptical mass distributions 
(SPEMDs) and for softened power-law elliptical potentials (SPEPs), and illustrate further the 
limitations of SPEPs. In §3 we present and simplify the quadrature solution of Schramm (1990) for 
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the deflection angle. We show that for the SPEMDs it is possible to approximate the integrand so 
that the integral can be done analytically. Although the result is a sum of series expansions, each 
series converges rapidly to high accuracy, even for mass densities with arbitrarily high ellipticity. 
We show how to similarly evaluate the magnification matrix. We also derive an expression for the 
gravitational potential, but it cannot be evaluated without a numerical integration. Finally, in §4 
we summarize our results. 



2. SPEMDs and SPEPs 



Consider a projected mass density £ at an angular diameter distance D d and a source at D s , 
with a distance of D^s from the lens to the source. Then the lens equation can be written as 



y = x - a(x) 



(1) 



where y is the source angle, x = {x\,x%) is the observed image angle, and a = {ot\,a.2) is the 
deflection angle scaled by a factor of Dd s /D s (see e.g. Schneider et al. 1992 for an introduction to 
gravitational lensing). The deflection is 



a(x) = W>(x), 



in terms of the potential 



^(x) = — f d 2 x' 'k(x') In |x — x'| 
vr J 

Equation |3] is the solution to the Poisson equation 

vV(x) = 2«(x) , 
where k = S/S cr in terms of the critical density 

c 2 D R 



The SPEMD is given by 



4vrG D d D ds 



p 2 + s 2 



E 2 



(2) 
(3) 

(4) 
(5) 

(6) 



where s is the core radius, rj is the power-law index, and E fixes the overall normalization. The 
dependence on position is through 



2 2,2/ 2 a 

p = X± + x 2 / cos p 



(7) 



where cos (5 is the axis ratio b/a and the ellipticity is e = 1 — cos f3. The parameter rj can vary from 
a modified Hubble profile (77 = 0) through isothermal (77 = 1) to a constant surface density sheet 
(rj = 2). 
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The SPEP is specified by a potential 



2El 



Pp ~l~ s p 



E 2 



where 



2 2,2/ 2 n 

p p = x 1 + x 2 / COS ftp 



The corresponding density as specified by the Poisson equation (f|) is 

x i + x \l cos4 Pp 



' pI + 4' 


'I 


E 2 
. p . 





pj + sj 



+ 1 + 



COS 2 /?„ 



(8) 



(9) 



(10) 



For p p 3> s p , the axis ratio of this density distribution is 

1 + (r/- l)cos 2 /3 p 



^4 P = cos /3p 



77 -sin /?„ 



i 

>7-2 



(11) 



(Equivalent forms for K p and A p have been derived by KK93 and Grogin and Narayan (1996), 
respectively). KK93 have shown that the isodensity contours acquire a concave part and become 
dumbbell shaped when cos 2 P v < 1 — r//3. Thus for a given r\ the requirement of convex density 
contours implies a smallest axis ratio A p that can be represented, or equivalently a maximum 
ellipticity of 



M = i 



2-2 



1 

V-2 



3 V 2, 

Figure 1 plots this maximum ellipticity as a, function of 77, cind. shows that 6 n 
0.5. An ellipticity greater than 0.65 cannot be modelled for any r\ < 2. 



(12) 

is typically around 



Assigning corresponding parameters between the SPEP and the SPEMD is somewhat arbitrary, 
but we follow KK93 and use the following conventions. For zero core radius and zero ellipticity, 
we require the SPEP and SPEMD to match, which yields the condition E = E p . We require the 
density axis ratios b/a far from the core to be equal, i.e. cos f3 = A p . Finally, we also require the 
central densities to match (if the core radii are nonzero), which yields 



1 

1-2 



1 + cos 2 Pp 

7] COS 2 P p 

Figure 1 shows that at a fixed rj, as the ellipticity is increased the contours of the SPEP become 
dumbbell shaped at the critical value e max (r/). This is illustrated in Figure 4 of KK93, which shows 
the density contours of several pairs of corresponding SPEMDs and SPEPs, for several values of 
b/a and with rj fixed at 1. Figure 2 illustrates the complementary case, where if we decrease rj at a 
fixed ellipticity the contours switch over from being convex to having the dumbbell shape. Figure 
2a shows the isodensity contours of an SPEMD with axis ratio b/a = 0.5. For simplicity, we set 
the core radius equal to zero. For the SPEMD, the contour shape depends only on the axis ratio 
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Fig. 1. — Maximum ellipticity representable by an SPEP model shown as a function of the power- 
law index rj. 



and is independent of rj. Figure 2b shows the isodensity contours of the corresponding SPEP with 
rj = 1.6 and the same axis ratio of 0.5 for the density contours. As we decrease rj keeping the axis 
ratio fixed, the SPEP contours remain convex until the value of rj = 1.25 shown in Figure 2c, but 
they become dumbbell shaped below this as illustrated in Figure 2d for rj = 0.4. 

Even with the problem of the SPEP regarding isodensity contour shape, it is of interest to 
directly compare the lensing behaviour of the two models, the SPEP and SPEMD. We compare 
them in the same configuration as KK93, i.e. for a source placed directly behind the center of the 
lens. In this case there are four images, with two images on the y-axis at (0, ±y) with magnification 
A y and two on the x-axis at (=tx, 0) with magnification A x . Figure 7 of KK93 shows the distance 
ratio x/y and the magnifications A x and A y for the SPEMD and corresponding SPEP, as a function 
of ellipticity for the singular (i.e. zero core) isothermal profile. In figure 3 we plot the distance ratio 
y/x and the magnification ratio A y /A x for singular profiles with rj = 0.5, rj = 1, and rj = 1.5. The 



-6- 



(a) b/a = 0.5 (b) b/a = 0.5 




SPEMD SPEP, 77= 1.6 



(c) b/a = 0.5 (d) b/a = 0.5 




SPEP, 77= 1.25 SPEP, 77 = 0.4 



Fig. 2. — Isodensity contours of SPEMDs and corresponding SPEPs. Panel (a) shows the SPEMD 
contours with axis ratio 0.5. For a fixed axis ratio, as we decrease 77 the SPEP contours switch 
from convex to dumbbell shaped, as shown in panels (b), (c), and (d). 



distance ratios behave similarly for the SPEMD (solid curves) and SPEP (dashed curves), though 
there are small differences at high ellipticity. The magnification ratio is an observable, not the 
individual image magnifications, and we find that the behavior of the magnification ratio at high 
ellipticity differs significantly between the SPEMD and SPEP even at r\ = 1 , but more so at smaller 
77. Thus while the SPEP produces some image configurations qualitatively similar to those of the 
SPEMD, quantitatively it is not an accurate substitute unless the ellipticity is small. 
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3. The deflection angle, magnification matrix, and potential of the SPEMD 

We begin with the solution of Schramm (1990) for any elliptical density of the form n{p), which 
we write concisely as 

ai(x 1 ,X2) = 2x\ cosp / -. — Tjdp , 

Jo xf + w 4 x| 

rp(xi,x 2 ) p' K (p')uj 3 

a 2 (xi,x 2 ) = 2x 2 cos/3 / — Tjdp , (13) 

Jo xf + w 4 X2 

A + r 2 + p /2 sin 2 /3 



A 2 = 



A + r 2 — p' 2 sin 2 /3 

!2 ■ 2 o , 2 2l 2 l/ i22 

p sm fj + X2 — xf + 4xf Xg 



where r 2 = x 2 + x 2 . The upper limit of integration p{x\,x 2 ) is given by equation (|7j). In what 
follows we can restrict x\ and X2 to be nonnegative without loss of generality, noting from equations 
( Pf ) that ai(xi,x 2 ) = ai(\x\\, \x 2 \) sign(xi) and a 2 (xi,x 2 ) = a 2 (\xi\, \x 2 \) sign(x 2 ). 

To simplify A, we switch variables to 

p' 2 sin 2 + x| - xf 

M = o • ( 14 ) 

2xix 2 



We find that u> 2 factorizes as lo 2 = (p + \J p 2 + 1) xl/x2, which is a key point for our method below 
to work. If we write k as a function of p, i.e. K,(p) = n(p), we obtain 

cos/VxiX2 -/ u 
a 1 {x l ,x 2 ) = r\- — / K(p)f(p)dp , 

sm p J [i 1 

cos/g^/xj^ /" M2 _, v., v, 

a 2 (xx,x 2 ) = t\- — / K(p)f(-p)dp , 

sm p J 



/( "> s \/7TT7?-?TT' (15> 



Ml 

f'2 



1 / X 2 Xi 

2 \xi x 2 , 

1 I x 2 x\ cos 2 /3 

2 I xi cos 2 (3 x 2 



This analysis so far is correct for any density of the form n(p). We now specialize to the 
SPEMD, which we can write as R{p) = q(p + s)~ 7 . In terms of our variables from §2, we have 
7 = 1 — 77/2, q = [2xix 2 /(E 2 sin 2 /3)]~ 7 , and s = —p\ + s 2 sin 2 (3/{2x\x 2 ). Defining v = p + s, we 
finally obtain 

ai(xi,x 2 ) = qh , 

a 2 (xi,x 2 ) = qh , (16) 
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in terms of a common factor 



Q = ^T~n Q ( 17 ) 

sin p 



and the integrals 

h = / v-y{v-s)du 

J 1/1 

/ 2 



/ V -if(s-v)du, (18) 



evaluated between the limits 



s 2 sin 2 /3 sin 2 /? /m x 2 , 
v\ = ~z 5 V 2 = v\ H — 1 t-j • (19) 



2xiX2 ' 2 \X2 XlCOS 2 /? 

To evaluate the integral I\ we make use of the fact that /(/x) depends only on the variable \x 
and we know its functional form, so we can approximate f(v — s) as a polynomial in either ^ or 1/u 
for various ranges of v from to oo (the lower limit v\ in ii is always nonnegative) . After some 
experimentation, we settle on a number of strategies for these polynomial expansions. 

For a small range of \i values around 0, e.g. ji a to we use Chebyshev polynomials (see e.g. 
Press et al. 1992) to construct a polynomial approximation for /(/i) of the form ^"=0^ c n^ n ■ Then 
we can evaluate the portion of I\ between v a = fi a + s and Vb = ^b + s to be (e.g. if 7 is not an 
integer) 



J Va 



n=N n=N m=n 



n=0 n=0 m=0 v 7 ' 



We can reduce the number of computations if s 3> \fib — fJ, a \, because we can then write = 
s _7 [l + {y— s)/s]~ 7 , expand this in powers of {y — s)/s, and use v — s as the variable of integration. 
The advantage of using Chebyshev polynomials is that the expansion converges rapidly if the range 
of /j, values is not too large. 

For /u > 1, we expand /(//) in inverse powers of fx, 

/( "» " 75^ - ■ <20> 

Consider just the first term (we deal similarly with the second term). We can write a contribution 
to I\ in this regime as 



I 



Vb dv 1 | -,- 7 -l/2 — ~ dU 



V2{v-sfn V2' ' Jv a /\s\ (P-n)3/2 ' 

where we have set v = v/\s\ and u = sign(s) is ±1. Then if v S> 1, we expand 

c-<.)-^ = ^(i + S + ^ + - 



- 9 - 



and then integrate. On the other hand, for small values of v up to a few we can again use Chebyshev 
polynomial approximation, this time for the function (v — u)~~ 3 / 2 , and then integrate (we do this 
separately for the two signs of u). In practice we find that it is better to perform an integration by 
parts and then use a polynomial approximation for the function {p — u) -1 / 2 instead. When u = +1 
we cannot extend this down to v = 1, so for a range of v near 1 we leave (p — 1)~ 3 / 2 unchanged 
but Taylor expand the term about v = 1. The result is an integrand which is a sum of powers 
of v — 1 and is easily integrated. When u = — 1, we find that we sometimes need high accuracy 
near v = 0, so for a small range of v <C 1 we Taylor expand (i> + l) -3 / 2 about the midpoint of this 
narrow range, and then integrate. 



We use the same set of methods when (— n) 3> 1, but here the expansion is 

V2 3 



-/i 4 v / 2(-Ai) 5/2 ' 



(21) 



In some cases we can use other methods to increase the speed. If the entire integration range 
of Ii is narrow, i.e. (1/2 — vi) "C 1, then we Taylor expand f{y — s) about the midpoint of the range 
and then integrate, without breaking up this range into the various regimes described above. 

The required accuracy determines how we specify the various regimes. We were motivated by 
the lens 0957+561 (e.g. Walsh et al. 1979, Young et al. 1981, Falco et al. 1991, Grogin & Narayan 
1996), where VLBI observations (Porcas et al. 1981, Gorenstein et al. 1988, Garret et al. 1994) 
require a relative accuracy in the deflection angle of one part in 10 5 . For algorithmic simplicity we 
prefer to use only two terms in the expansions ( f20| ) and (^Tj) , so this high accuracy forces us to use 
them only for > 15.7. For —15.7 < \i < 15.7 we use the direct Chebyshev expansion of /(//). 
We find it best to break up this range into 7 intervals and construct a Chebyshev expansion in 
each, with the number of terms in the expansions ranging from 6 to 9, keeping the relative accuracy 
roughly constant. 

The integral I2 is evaluated using exactly the same expansions as I\ but with /(— /u) substituted 
for /(/i). The components of the 2x2 inverse magnification matrix <5jj — doti/dxj also yield similar 
integrals. Using equations ( |l6|) through (|l9| ) we derive 



da.\ 
dx 2 



V 1 



Oil ~ 

— + q 



-7,, ^ dv\ j ds 
1 dxi 3 dx\ 



1 

2 " 
dai 
dx 2 



7 



+ q 



a 2 

x 2 
1 

2" 7 



v 2 7 /(s - v 2 ) 



dv 2 

dx 2 



7 f(s 



x 2 



-\ du2 

f{V2 - S) dx~ 2 



. dv\ ds 
ox 2 ax 2 

f^- s )l^r-hT— 



(22) 



dx 



dx 2 



Just as with I\ and I 2 , we evaluate the integrals 



h 
h 



b>2 



v 1 f ' {y — s)dv , 
z/ _7 /'(s — v^dv , 



(23) 
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where we approximate f'(p) = df(n)/dfi in the same ways as we created polynomial approximations 
for /(//) above. The other quantities appearing in equations (22) can be evaluated (for i = 1 and 
2) as 



dxi 
dxi 



v\ ds 

x 2 dx 2 

dv\ sin 2 (5 \x\ 

dxi 2xi [x 2 



dv x 1 
dx^ —j i 

X2 



x\ cos 2 (3 



Xi X2 
X2 Xi 



-iy 



(24) 



If we apply the same amount of computational effort to I3 and I4 as for 1% and I2, we never- 
theless find that the relative accuracy is lower for the magnification det _1 |5jj — dcti/dxj\ than it 
is for the components of a, for a couple reasons. First, the various polynomial expansions tend to 
converge more slowly for f'(p) than for f{p); and second, the components of a are directly propor- 
tional to l\ and I2, but in equations (^) we sometimes have to subtract nearly equal terms and 
get a low accuracy for doti/dxj even if ^3 and I4 are evaluated accurately. In practice we achieve a 
maximum relative error of 5 x 10~ 6 for a and 6 x 10~ 4 for the magnification, and a typical relative 
error of 1 x 10~ 6 for a and 5 x 10~ 5 for the magnification. A higher accuracy is not needed for the 
magnification, given typical measurement errors on fluxes. 

In terms of running time, with the above accuracy we can evaluate the deflection angle and 
magnification matrix of the SPEMD roughly 20 times faster than the brute force method (which 
requires 5 integrals), although this is still about 15 times slower than the SPEP, or any other model 
with the deflection given by a simple formula. This speedup should make the SPEMD model useful 
for applications in which repeated numerical integrations made it previously unusable. 

To evaluate the gravitational part of the time delay, we also need to evaluate the potential 
?/>(x). Schramm (1990) gives a quadrature expression for the potential of a single elliptical shell. 
We obtain the potential of any elliptical density n{p) by integrating over shells, 

rp(x u x 2 ) n(p') du 

= cos/3 J p/=Q y/{p/ 2 + u){pl 2 C0S 2 (3 + u) P <P ) d P ■ ( 25 ) 

where 

7(p') = ^p-p ,2 (l + cos 2 /?) + A} (26) 
in terms of the last equation of (|l3|). Performing the inner integral, we derive 

tp(xi,X2) = 2cos/3 / In 

Jp'=0 

In this case the integrand does not factorize when we substitute equation (0), so the potential 
must be numerically integrated and cannot be speeded up. In most applications, though, it is not 
necessary to evaluate the potential a large number of times. 



A + r 2 - p' 2 sin 2 + \J A + r 2 + p 12 sin 2 (5 
v^p'Cl + cos 0) 



p'nip'W . (27) 



4. Conclusions 



A mass density profile with elliptical isodensity contours is a natural lens model to try when 
axisymmetric models fail. Previously, however, the easiest way of evaluating the deflection has been 
to numerically integrate the solutions of Schramm (1990). After simplifying these solutions we have 
shown that for the family of SPEMD mass distributions, the deflection angle and magnification 
matrix can be evaluated very fast and with high accuracy. Our implementation achieves a relative 
accuracy of 5 x 10~ 6 in the deflection and 6 x 10~ 4 in the magnification while running 20 times 
faster than a procedure based on the numerical integrations. We have also derived an expression 
for the potential, although this quantity must be numerically integrated. 

As noted by Schneider & Weiss (1991), combinations of two or more SPEMDs with different 
parameters can be used to construct more general density profiles with several scales. We thus 
expect SPEMDs to be more widely used, particularly for cases of high ellipticity in which the 
alternative SPEPs develop the artificial feature of dumbbell-shaped contours. 

I thank Ed Turner for valuable discussions. This work was supported by Institute Funds. 
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Fig. 3. — Properties of the four images of a source placed directly behind the lens center, for 
SPEMDs and corresponding SPEPs. Several different powerdaw indices rj are considered. The 
panels on the left show the distance y of the images on the y-axis over the distance x of the images 
on the x-axis, and the panels on the right show the magnification ratio A y /A x of the two types 
of images. In each case the appropriate quantity is plotted as a function of ellipticity, with the 
SPEMD as the solid curve and the SPEP as the dashed curve. 



